2D or 3D FDTD Modeling of Photonic Crystal Waveguides?

The aim of this paper is to show how 2D approaches, often used to investigate numerically the real devices as they are faster and simpler respect to 3D techniques, may lead to wrong results in modeling Photonic Crystal structures. As a test case, the determination of the so called Mini Stop Band (MSB) for the fundamental mode in a W1 waveguide is considered. A MSB is found when the structure, reduced to a 2D one, is analyzed by either Finite Difference Time Domain (FDTD) or Plane Wave Method approaches, but disappears when the real device is investigated using the FDTD-3D technique. Moreover, discrepancies exist between the 2D and the 3D transmission curves. The 3D approach is then necessary to get accurate results and determine the structure correct behavior. Introduction Photonic Crystals (PC) are periodic arrangements of dielectric or metallic materials. In these years they have been intensively studied to design and realize many optical components, exploiting their ability to control light propagation [1]. A single defect in a PC can localize the electromagnetic field in the frequency Band Gap (BG) of the unperturbed structure [2], while a sequence of defects can guide an optical signal along it, creating frequency-selective waveguides. Many numerical approaches have been proposed so far to study these phenomena either in the frequency or in the time domain [1, 3]. Among them, the FDTD [3] (Finite Difference Time Domain) techniques and the Plane Wave Method (PWM) are probably the most popular. Their advantages and problems have been evidenced, for example, comparing results obtained with the PWM and the FDTD in a 2D case for a single defect [1] or with FDTD-2D and FDTD-3D for a line waveguide [4]. Generally speaking, 2D approaches are attractive as they are fast and not numerically intensive. However, they are not expected to guarantee the same accuracy of 3D techniques. In this paper we will show that, in effect, though 2D modeling is useful for qualitative analysis, for example to calculate features of newly conceived devices, 3D models cannot be avoided for quantitative device design. To this purpose, we will compare, in the next sections, results obtained using PWM, FDTD-2D and FDTD-3D on a test structure based on a PC with triangular lattice of air holes carved in a 3 layer dielectric slab. The 3D structure will first be reduced to a 2D one by the effective index method (EIM) approximation [5]. The fundamental mode of the 2D device, studied using both FDTD and PWM [6], will be shown to present a certain number of so called Mini Stop Bands (MSB) [6,7], some of which, however, disappear using a FDTD-3D analysis. Differences appear also in the position of minima of the 2D transmission coefficient. This will lead to conclude that 3D simulations are necessary at least in the final device design. Device characteristics and numerical methods The studied PC has a triangular lattice with lattice constant a = 0.215 μm and air holes with radius = 0.1032 μm. The PC is carved in a three layer slab with values of the relative dielectric constants equal to 11.5, 13.6 and 11.5. This structure exhibits a gap for modes with electrical field parallel to the cylinder axes (TM modes). As anticipated, the real 3D structure can be reduced to a 2D one by the EIM approximation. This is equivalent to substitute the layered structure by a homogeneous substrate with refractive index equal to the effective index of the slab. We then analyzed using both PWM and FDTD-2D: (i) the BG for the perfect PC, (ii) the dispersion curves and (iii) the pattern of the field of a so called W1 waveguide (which can be obtained imply removing one row of holes). The FDTD-2D we used is based on the Yee’s algorithm [8] and can calculate the transmission coefficient and the dispersion curves using the supercell method [9]. The geometry of the studied structure is plotted in fig. 1. The computational domain consists in 372x623 cells with steps ∆x = ∆y = 0.005 μm and includes 8 GT-PML layers for each boundary. Progress in Electromagnetic Research Symposium 2004, Pisa, Italy, March 28 – 31 382 Figure 1. Supercell for the FDTD-2D calculation. Parameters used in the simulation are radius/a=0.48, air holes in dielectric medium with effective index= 3.61 The spatial distribution of the excitation field is a summation of plane waves [9], pulsed in the time domain, satisfying the Bloch conditions in each primitive cell. As boundary conditions, the sine-cosine method [3] is implemented for each wall of the supercell. The time evolution of the field is sampled in several waveguide points to avoid missing of some mode contributions, Fourier transformed and results averaged. The peaks in the spectrum represent the eigenmodes of the structure. Dispersion curves are obtained plotting the frequency values of the obtained peaks versus the wave vectors of interest for the waveguide. The used computational domain consists in 25×260 cells. The FDTD-2D results have been compared with those obtained via PWM using a 1×8 supercell with 2401 plane waves. The real 3D structure, on the contrary, has been studied using a parallel version of the FDTD3D simulator running on a cluster of personal computers [10] realized to overcome problems related to longer CPU time and larger memory requirements of the 3D model. In the 3D case, the transmission coefficient of the fundamental mode of the W1 waveguide (made in the three layer slab) has been calculated exciting the structure with a ridge waveguide, whose mode was independently determined by a Finite Difference modal solver. The mesh used in the 3D simulations consists in 372×258×100 nodes, with ∆x = ∆y = 0.005 μm, ∆z = 0.008 and 8 GT-PML layers for each boundaries. To perform a wide band analysis, and obtain the frequency response, the input field is temporally modulated by a Gaussian pulse with spectrum covering the frequency band of interest. Then, the transmission coefficient is calculated comparing, at each frequency, the Poynting vectors at the output and input sections of the PC waveguide. In particular, the incident and transmitted electromagnetic fields on these two sections are Fourier transformed and used to compute, for each harmonic, the flux of the Poynting vector in each single cell. The total incident and transmitted powers are then computed averaging, on the whole extension of the I/O sections, those “local” values.